﻿using System;
using System.Collections.Generic;
using System.Text;
using System.IO;

namespace GeoFly
{
    //土壤参数结构
    public class SoilInfo
    {
        public SoilInfo()
        {            
            Horton_K = 0.05;
            albedo = 0.23;
            
        }
        /// <summary>
        /// 创建土壤层结构
        /// </summary>
        /// <param name="LayerCount">土壤层数</param>
        public void CreateLayers(int LayerCount)
        {
            this.SoilLayers = new SoilLayer[LayerCount];
            this.LayerCount = LayerCount;
            for (int i = 0; i < LayerCount; i++)
            {
                SoilLayers[i] = new SoilLayer();
            }
        }
        /// <summary>
        /// 土壤结构初始化和输入、输出处理
        ///从文件中加载每一种指定土壤类型的主要输入物理参数（固有参数）
        /// </summary>
        /// <param name="FileName"></param>
        /// <returns></returns>
        public void ReadSoilFile(string FileName)
        {
            StreamReader sr = new StreamReader(FileName);
            string[] strs= sr.ReadLine().Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
            this.Soil_Name = strs[1];

            strs=sr.ReadLine().Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
            this.LayerCount = Convert.ToInt16(strs[1]);

            strs = sr.ReadLine().Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
            this.rootdepth = Convert.ToDouble(strs[1]);

            strs = sr.ReadLine().Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
            this.albedo = Convert.ToDouble(strs[1]);

            strs=sr.ReadLine().Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
            this.Horton_K = Convert.ToDouble(strs[1]);

            strs=sr.ReadLine().Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
            this.InitSWP = Convert.ToDouble(strs[1]);

            sr.ReadLine();
            sr.ReadLine();

            this.SoilLayers = new SoilLayer[LayerCount];
            for (int i = 0; i < LayerCount; i++)
            {
                this.SoilLayers[i] = new SoilLayer();
                strs=sr.ReadLine().Split(new char[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
                this.SoilLayers[i].SL_Z = Convert.ToDouble(strs[2]);
                this.SoilLayers[i].SL_BD = Convert.ToDouble(strs[3]);
                this.SoilLayers[i].SL_AWC = Convert.ToDouble(strs[4]);
                this.SoilLayers[i].SL_Sat_K = Convert.ToDouble(strs[5]);
                this.SoilLayers[i].SL_Stable_F = Convert.ToDouble(strs[6]);
                this.SoilLayers[i].SL_Org_C = Convert.ToDouble(strs[7]);
                this.SoilLayers[i].SL_Clay = Convert.ToDouble(strs[8]);
                this.SoilLayers[i].SL_Silt = Convert.ToDouble(strs[9]);
                this.SoilLayers[i].SL_Sand = Convert.ToDouble(strs[10]);
                this.SoilLayers[i].SL_Rock = Convert.ToDouble(strs[11]);
                this.SoilLayers[i].SL_Init_F = Convert.ToDouble(strs[12]);
            }
            this.CalcSoilPara();
        }

        public bool WriteSoilFile(string FileName)
        {
            throw new Exception("");
        }	//将土壤类型的主要物理参数保存成土壤参数文件

        //土壤参数计算
        private void CalcSoilPara()
        {
            for (int i = 0; i < LayerCount; i++)
            {
                this.SoilLayers[i].CalcSoilPara();
            }

            double acumudepth, sumpor, lyrdepth, pormm, stainfil, stainit, suminfil, suminit, sumsat_k, sat_k;
            acumudepth = sumpor = lyrdepth = pormm = stainfil = stainit = suminfil = suminit = sat_k = sumsat_k = 0;
            
            for (int i = 0; i < LayerCount; i++)
            {
                lyrdepth =SoilLayers[i].SL_Z - acumudepth;
                pormm =SoilLayers[i].SL_Por * lyrdepth;
                sat_k =SoilLayers[i].SL_Sat_K * lyrdepth;
                stainfil =SoilLayers[i].SL_StaInfil * lyrdepth;
                stainit =SoilLayers[i].SL_InitInfil * lyrdepth;

                sumpor += pormm;
                suminfil += stainfil;
                suminit += stainit;
                sumsat_k += sat_k;

                SoilLayers[i].SL_Sat = (SoilLayers[i].SL_Por - SoilLayers[i].SL_WpRatio) * lyrdepth;
                this.SP_Sat += SoilLayers[i].SL_Sat;
                SoilLayers[i].SL_Wp = SoilLayers[i].SL_WpRatio * lyrdepth;                
                SP_Wp +=SoilLayers[i].SL_Wp;
                SoilLayers[i].SL_P_Fc = (SoilLayers[i].SL_FcRatio -SoilLayers[i].SL_WpRatio) * lyrdepth;
                if (SoilLayers[i].SL_P_Fc < 0)
                    throw new Exception("参数错误");
                SP_Fc +=SoilLayers[i].SL_P_Fc;
                SoilLayers[i].SL_SW =SoilLayers[i].SL_P_Fc * InitSWP;
                SoilLayers[i].SL_HK = (SoilLayers[i].SL_Sat -SoilLayers[i].SL_P_Fc) /SoilLayers[i].SL_Sat_K;
                if (SoilLayers[i].SL_HK < 1)
                    SoilLayers[i].SL_HK = 1;
                Initial_SP_Sw +=SoilLayers[i].SL_SW;
                acumudepth =SoilLayers[i].SL_Z;
            }
    
            SP_Por = sumpor /SoilLayers[LayerCount-1].SL_Z;
            SP_BD = 2.65 * (1 - SP_Por);
            SP_Stable_Fc = suminfil / SoilLayers[LayerCount - 1].SL_Z;
            SP_Init_F0 = suminit / SoilLayers[LayerCount - 1].SL_Z;
            SP_Sat_K = sumsat_k / SoilLayers[LayerCount - 1].SL_Z;
            if (Initial_SP_Sw < SP_Wp)
                Initial_SP_Sw = SP_Wp * 1.05;										//	FRZZ 2006-02-06 

            double f = 6.5309 - 7.32561 *SoilLayers[0].SL_Por + 0.001583 * Math.Pow(SoilLayers[0].SL_Clay, 2)
                + 3.809479 * Math.Pow(SoilLayers[0].SL_Por, 2) + 0.000344 *SoilLayers[0].SL_Sand *SoilLayers[0].SL_Clay
                - 0.049837 *SoilLayers[0].SL_Por *SoilLayers[0].SL_Sand
                + 0.001608 * Math.Pow(SoilLayers[0].SL_Por, 2) * Math.Pow(SoilLayers[0].SL_Sand, 2)
                + 0.001602 * Math.Pow(SoilLayers[0].SL_Por, 2) * Math.Pow(SoilLayers[0].SL_Clay, 2)
                - 0.0000136 * Math.Pow(SoilLayers[0].SL_Sand, 2) * SoilLayers[0].SL_Clay
                - 0.003479 * Math.Pow(SoilLayers[0].SL_Clay, 2) *SoilLayers[0].SL_Por
                - 0.000799 * Math.Pow(SoilLayers[0].SL_Sand, 2) *SoilLayers[0].SL_Por;

            SP_WFCS = 10 * Math.Exp(f);
            this.TPercolation = (this.SP_Sat - this.SP_Fc) / this.SP_Sat_K;
        }
        /// <summary>
        /// 获取土壤的初始含水量，在执行该属性前必须先执行过CalcSoilPara()方法
        /// </summary>
        public double Initial_SP_Sw=0;
        //利用固有土壤参数计算其它有用的土壤物理参数
        /// <summary>
        /// 计算干旱指数	
        /// </summary>
        /// <param name="bConst"></param>
        public double CalcAridIndex(bool bConst,double SP_Sw) 
        {
            if (bConst)		//用常数计算
                return (SP_Fc - SP_Wp) / SP_Fc;
            else			//用变量计算
                return  (SP_Sw - SP_Wp) / SP_Fc;
           
        }
        public double Initial_Arid
        {
            get
            {
                if(Initial_SP_Sw<=0 || SP_Wp<=0 || SP_Fc<=0)
                    throw new Exception("土壤参数可能没有初始化");
                return (this.Initial_SP_Sw - SP_Wp) / SP_Fc;
            }
        }
        /// <summary>
        /// 土壤水赤字比例	
        /// </summary>
        /// <returns></returns>
        public double SoilWaterDeficitPercent(double SP_Sw) 
        {
            return (1.0 - SP_Sw / SP_Fc) * SP_Por;
        }
       /// <summary>
       /// 求土壤水赤字含水量
       /// </summary>
       /// <param name="SP_Sw">实际土壤含水量</param>
       /// <returns></returns>
        public double SoilWaterDeficitContent(double SP_Sw) 
        {
            if (SP_Sw <= 0 || SP_Fc <= 0 || SP_Por <= 0 || rootdepth <= 0)
                throw new Exception("数据不能为负");
            return (1.0 - SP_Sw / SP_Fc) * SP_Por * rootdepth;
        }
        /// <summary>
        /// 土壤剖面平均含水量
        /// </summary>
        /// <returns></returns>
        public double SoilAvgWater(double SP_Sw)
        {
            if(SP_Sw<=0 || rootdepth<=0)
                throw new Exception("数据不能为负");
            return SP_Sw / rootdepth * 100;
        }


        //外部输入参数
        /// <summary>
        /// 土壤名称
        /// </summary>
        public string Soil_Name;
        public int Soil_id;
        /// <summary>
        /// 土壤分层数
        /// </summary>
        public int LayerCount;
        /// <summary>
        /// 植被根深(mm)
        /// </summary>
        public double rootdepth=0;

        public double albedo=0;
        public double Horton_K=0;
        public double InitSWP=0;

        public SoilLayer[] SoilLayers;
        
        /// <summary>
        /// 土壤剖面平均湿容重
        /// </summary>
        public double SP_BD=0;
        /// <summary>
        /// 土壤剖面平均孔隙率
        /// </summary>
        public double SP_Por=0;
        /// <summary>
        /// 田间持水量时土壤含水量
        /// </summary>
        public double SP_Fc=0;
        /// <summary>
        /// 土壤饱和时土壤含水量
        /// </summary>
        public double SP_Sat=0;
        /// <summary>
        /// 枯萎点时土壤含水量
        /// </summary>
        public double SP_Wp=0;
        /// <summary>
        /// 初始土壤实际含水量
        /// </summary>
        //public double SP_Sw;
        /// <summary>
        /// 
        /// </summary>
        public double SP_WFCS=0;
        /// <summary>
        /// 土壤剖面稳定下渗率
        /// </summary>
        public double SP_Stable_Fc=0;
        public double SP_Init_F0=0;
        //public double Initial_Arid=-1;
        /// <summary>
        /// 土壤剖面平均饱和导水率
        /// </summary>
        public double SP_Sat_K=0;
        public double TPercolation=0;
    };
}